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Suitable gauge conditions are fundamental for stable and accurate numerical-relativity simulations of inspi- 
ralling compact binaries. A number of well-studied conditions have been developed over the last decade for 
both the lapse and the shift and these have been successfully used both in vacuum and non- vacuum spacetimes 
when simulating binaries with comparable masses. At the same time, recent evidence has emerged that the 
standard "Gamma-driver" shift condition requires a careful and non-trivial tuning of its parameters to ensure 
long-term stable evolutions of unequal-mass binaries. We present a novel gauge condition in which the damp- 
ing constant is promoted to be a dynamical variable and the solution of an evolution equation. We show that this 
choice removes the need for special tuning and provides a shift damping term which is free of instabilities in 
our simulations and dynamically adapts to the individual positions and masses of the binary black-hole system. 
Our gauge condition also reduces the variations in the coordinate size of the apparent horizon of the larger black 
hole and could therefore be useful when simulating binaries with very small mass ratios. 

PACS numbers: 04.40.-b,04.40.Dg,95.35.+d 



I. INTRODUCTION 

Five years after the first demonstrations iflJO that the nu- 
merical solution of the inspiral and merger of binary black 
holes (BBHs) was within the technical and computational 
capabilities of many numerical-relativity groups, our under- 
standing of this process has expanded beyond the most opti- 
mistic predictions. Numerical-relativity simulations of black- 
hole binaries have been performed in a large region of the 
possible space of parameters (see [4, 5] for two recent re- 
views). In addition, the results of these simulations have been 
exploited on several fronts. In gravitational wave data anal- 
ysis, they have been used to produce template banks that in- 
crease the distance reach of detectors [6-8] and to aid the cal- 
ibration of search pipelines ll9 fflTI . In astrophysics they have 
been used to determine the properties of the final black hole 
(BH) of a BBH inspiral and merger (see lfT2UT4l and refer- 
ences therein for some recent work) and hence assess the role 
that the merger of supermassive BHs plays in the formation of 
galactic structures [ 15 , 16]. In cosmology they have been used 
to study the electromagnetic counterparts to the merger of su- 
permassive BH binaries and hence deduce their redshift IT71 - 

sa. 

Despite this extensive and rapid progress, there are por- 
tions of the space of the parameters that still pose challenges 
for numerical simulations, in particular those involving bi- 
naries with BHs that are maximally spinning (but see ETI 
for some progress in this direction) or with small mass ra- 
tios q = M2/A/1 < 1. This latter problem is potentially a 
rather serious one since the computational costs scale in gen- 
eral quadratically with the inverse of the mass ratio of the sys- 
tem. The inspiral timescale is inversely proportional to q (the 
smaller BH spends more orbits per frequency interval during 
its inspiral onto the larger BH) and the timestep limit in ex- 
plicit numerical schemes also decreases linearly with the BH 
size and hence inversely proportionally with q. 

While some progress has been made recently when simu- 
lating binaries with mass ratios as small as q — 1/10 l22ll23l . 



it is clear that some significant technical changes are needed 
in order to tackle mass-ratios which are much smaller. One 
of these changes may consist in adopting implicit numeri- 
cal schemes in which the timestep limitation is set uniquely 
by the truncation error and not by the smallest spacing of 
the spatial numerical grid. Another improvement could come 
from the use of better spatial gauge conditions that, by bet- 
ter adapting the coordinates to the different curvatures of the 
spatial slice, may reduce the numerical error and hence the 
computational cost. Recent numerical simulations have re- 
vealed that the standard "Gamma-driver" shift condition E4ll 
requires a careful tuning of its parameters to ensure long-term 
stable BBH evolutions when considering unequal-mass bina- 
ries. Work to alleviate some of these problems has been re- 
cently started fl25H27jl and has so far concentrated on adapting 
the damping term 7/ in the shift condition (see Sect.|II]for a de- 
tailed discussion of the gauge and of the damping term) to bet- 
ter suit the uneven distribution of local curvature on the spatial 
hypersurface as the binary evolves. In practice, while a con- 
stant value of 77 has worked well for comparable-mass bina- 
ries, the investigations reported in [25-27 1 have suggested the 
use of damping factors that have a spatial dependence adapted 
to the location of the BHs. The first non-constant prescrip- 
tion for 77 l25l used an expression which adapts to the mass 
of each BH via the conformal factor, but it was found to lead 
to large errors at mesh refinement boundaries [26 1 and so was 
replaced by a simpler analytical form containing constant pa- 
rameters which need to be tuned to the masses of the BHs. 
At present it is not clear whether the tuning made with mass 
ratios q > 0.25 will be effective also for much smaller mass 
ratios. 

In this paper we propose a different approach to the prob- 
lem of a dynamical damping term for symmetry seeking shift 
conditions and promote 7/ to be a fully dynamical variable, as 
has proven very effective for the shift vector and for the T 
variables. The resulting gauge condition dynamically adapts 
to the individual positions and masses of the BBH system 
without the need for special tuning and remains well-behaved 
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(i.e., smooth and bounded) at all times. Considering three dif- 
ferent options for the source term in the evolution equation 
for r\, we show that an expression which is very simple to 
implement leads to numerical errors which are comparable to 
or smaller than the constant rj case. Furthermore, this choice 
reduces the dynamics in the coordinate size of the apparent 
horizon of the smaller BH and could therefore be useful when 
simulating binaries with very small mass ratios. 

The structure of the paper is as follows. In Sect. [TT] we 
summarise the numerical infrastructure and the mathemati- 



cal setup used in our simulations, while Sect. Ill is dedi- 



cated to a review of the slicing and spatial gauge conditions 
and to the discussion of our novel approach. Sections IV 
and [V] are dedicated to the discussion of the results of apply- 
ing the new gauge to simulations of single nonspinning BHs 



(Sect. IV i and to systems with BHs having either equal or 
unequal masses (Sect. [Vj. Finally, the conclusions and the 
prospects for future work are detailed in Sect. VI We use 
a spacelike signature (— , +, +, +) and a system of units in 
which c = G = M Q = 1. 




II. NUMERICAL SETUP 



FIG. 1: Schematic diagram of a typical Llama grid setup in the 
(a;, y) plane. Note the inner Cartesian grid with box-in-box AMR 
which is joined to a 6-patch multiblock structure (only 4 of these 
patches are shown in the (x, y) plane). 



The numerical setup used in the simulations presented here 
is the same one discussed in 11281 and more recently applied 
to the Llama code described in (29). The latter makes use 
of higher-order finite-difference algorithms (up to 8th order 
in space) and a multi-block structure for the outer computa- 
tional domain, which allows one to move the outer boundary 
to a radius where it is causally disconnected from the binary. 
We refer the reader to the papers above for details and here 
simply note that we solve the Einstein equations in vacuum 
with a conformal and traceless formulation of the equations, 
in which the conformal factor has been redefined as W = 
[(det(7 a f,)] -1 / 6 , or in terms of the metric j a b = W 2 ~f a b- The 
corresponding evolution equation is therefore 

d t W - @%W = ^WaK - ^W&F , (1) 

where K is the trace of the extrinsic curvature (see 11281 
and [29 1 for details on our specific implementation). 

The computational infrastructure of the Llama code is 
based on the Cactus framework [30, 31] and the Carpet 
ll32l |33l mesh-refinement driver, and implements a system 
of multiple grid patches with data exchanged via interpola- 
tion E9l . We use a central cubical Cartesian patch containing 
multiple levels of adaptive mesh refinement, with higher res- 
olution boxes tracking the location of each BH, i.e., "moving 
boxes". This is surrounded by 6 additional patches with the 
grid points arranged in a spherical-type geometry, with con- 
stant angular resolution to best match the resolution require- 
ments of radially outgoing waves. This allows us to evolve 
out to very large radii at a tiny fraction of the computational 
cost which would be necessary to achieve the same resolution 
with a purely Cartesian code. We use one patch for each of the 
±x, ±y and ±z axes, which leads to an inner spherical inter- 
patch boundary and to a spherical outer boundary. The latter 



can be placed at very large distances and we choose it to be it 
causally disconnected from the surfaces on which we compute 
the waveforms for the duration of the simulation. Figure [T] 
shows a schematic diagram of a typical Llama grid setup in 
the (x, y) plane. Note the inner Cartesian grid with moving- 
boxes which is joined to a 6-patch multiblock structure (only 
4 of these patches are shown in the (x, y) plane). The detailed 
grid structure used in each run is listed in Table [IJ and for the 
purpose of comparison, the resolution of each simulation is 
indicated by the grid spacing ho /M of the coarsest Cartesian 
grid. The unit M is chosen such that each BH has mass 0.5M 
in both the single and binary BH cases. In all cases the coars- 
est resolution is also equal to the radial spacing in the angular 
patches. 

III. GAUGE CONDITIONS AND DYNAMICAL DAMPING 
TERM 

As mentioned in the introduction, for the formulation of the 
Einstein equations we adopt, the use of suitable gauge con- 
ditions was the last obstacle to overcome in order to obtain 
long-term stable simulations of BBHs O [3] . In what is now 
the standard moving-puncture recipe, the lapse a is evolved 
using a singularity-avoiding slicing condition from the 1 + log 
family l34l 

d t a-j3 i d i a = -2aK, (2) 

while the shift vector j3 l is evolved using the hyperbolic 
Gamma-driver condition [24| 

d t /3 a - p%p = \B a , (3) 
d t B a -fi l d l B a = d t f a - /3 l dif a ~f]B a . (4) 
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Configuration 


h /M 


-/V ang . 


Rin/M Rout/M JViev. 


n/M 


single BH 


0.96 


21 


39.36 400.00 6 


(12,6,3,1.5,0.6) 


BBH, 3 = 1 


0.96 


21 


39.36 1980.48 6 


(12,6,3,1.5,0.6) 


BBH, g = l/4 


(0.8,0.96,1.12) 


(23,27, 33) 


49.92 2545.92 (6,8) 


(12,6,3,1.5,0.8,0.4,0.2) 



TABLE I: Numerical grid parameters of each BH configuration studied, ho is the grid spacing on the coarsest Cartesian grid, which is equal 
in all cases to the radial grid spacing in the angular patches. iV ang . is the number of cells in the angular directions in the angular patches. R^ n 
and i?out are the inner and outer radii of the angular patches. JVi ov . is the number of refinement levels (including the coarsest) on the Cartesian 
grid, and r\ indicates that a cubical refinement box of side 2r\ is centred on the BH on level "1", with level being the coarsest (in the case that 
the boxes overlap, they are replaced with a single box enclosing the two). The unit M is chosen such that each BH has mass 0.5A/ in both the 
single and binary BH cases. 



We recall that the Gamma-driver shift condition is similar to 
the Gamma-freezing condition d t T k = which, in turn, is 
closely related to the minimal distortion shift condition (35). 
The differences between these two conditions involve the 
Christoffel symbols and are basically due to the fact that the 
minimal distortion condition is covariant, while the Gamma- 
freezing condition is not (see the discussion in [36 1). 

The coefficient rj of the last term in Q is usually referred 
to as the damping term and plays a fundamental role in our 
investigation. It was originally introduced to avoid strong os- 
cillations in the shift and experience has shown that by tuning 
its value, it is possible to essentially "freeze" the evolution of 
the system at late times f37l . In simulations of inspiralling 
compact binaries this damping term is typically set to be con- 
stant in space and time and equal to 2 /M for BBHs and equal 
to 1/M for binaries of neutron stars, where M is the sum 
of the masses of the BHs or neutron stars [38 39]. Similar 
values have also been shown to yield stable evolutions in the 
case of mixed binaries with mass ratios q ~ 1/6 [40]. While 
this choice works well for binaries with comparable masses, a 
simple dimensional argument shows that it will cease to be a 
good one for binaries with unequal masses. Since the dimen- 
sion of rj is inverse mass, and the relevant mass is the mass of 
each individual BH, as the BH mass decreases a larger value 
of rj will be needed to maintain a similar damping effect. For 
an unequal mass system, this is impossible with a constant rj. 



1. Position-Dependent Damping Term 

To overcome the limitations imposed by a constant-in- 
space damping term, various recipes have been proposed re- 
cently in the literature. A first suggestion was presented in 
Ref. Il25l . where the damping term was specified by the func- 
tion 



y/^diWdjW 



(5) 



where j v is the inverse of the conformal 3-metric and R is 
a dimensionless constant chosen in such a way that RqAIbti 
corresponds to the Schwarzschild radial coordinate for the sta- 
tionary state, i.e., Rq ~ 1.31241 flTll . However, as discussed 
subsequently by the same authors [26], expression |5]l leads 
to sharp features (spikes) which can produce coordinate drifts 



and affect the stability of the simulations. To remove these 
drawbacks, alternative forms were suggested that read respec- 
tively 1 26 1 



Vmgb = A + 



l + wi(f?)« 1 + w 2 (ff) r 



(6) 



and 



?7mgb 



= A + Cie- Wli ^" +C 2 e- W2i ^ )n , (7) 



where w± and W2 are positive parameters chosen to change 
the width of the functions based on the masses of the two 
BHs. The power n is a positive integer which determines the 
fall-off rate, while the constants A, C\, and C<i are chosen 
to provide the desired values of rj at the punctures and at in- 
finity. Finally, the dimensionless radii f\ and fi are defined 
as fi = \fi — r]/|r*i — r%|, where i is either 1 or 2, and fi is 
the position of the i-th BH. In addition to the new prescrip- 
tions |6|-(|7]), which provide appropriate values both near the 
individual punctures and far away from them with a smooth 
transition in between, there is evidence that they also lead to 
a smaller truncation error. In particular, when examined for 
toi = W2 = 12, n = 1, the waveforms produced when us- 
ing Eq. |6]l showed less deviation with increasing resolution 
than using a constant rj. Similar results were found when us- 
ing Eq. (FT), leading the authors to the conclusion that the pre- 
ferred definition for the damping term is Eq. (|6]l because it is 
computationally less expensive. 

It was shown recently E71 for the Gamma-driver shift con- 
dition that there is a stability limit on the time step size which 
depends on rj. This limitation comes from the time integrator 
only and is not dependent on spatial resolution. One of the 
proposed solutions to this problem is to taper rj with a func- 
tional dependence of the type 



R 2 



R? 



(8) 



where r is the coordinate distance from the centre of the BH 
and R is the radius at which one makes the transition between 
an inner region where 775 is approximately equal to ?y and an 
outer region where rjs gradually decreases to zero. We find 
that this form does indeed help in removing potential insta- 
bilities and, as we will discuss in the following section, can 
be used also for a fully dynamical definition for the damping 
term. 
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2. Evolved Damping Term 

While the analytic prescriptions discussed in the previous 
section have been shown to be effective when suitably tuned, 
it is not clear whether they will be equally effective for dif- 
ferent mass ratios, nor how to choose the parameters without 
case-by-case tuning. In view of this and in order to derive 
an expression which adapts dynamically in space and time to 
the local variations of the shift vector, we have promoted the 
damping term to be an evolved variable with the simple equa- 
tion 



d t i] - P l dir] 



1 

M 



(-V + S(r)) 



(9) 



The function S(r) is a position-dependent source term which 
can be chosen freely. The first term on the right-hand-side 
is introduced to induce an exponential decay of the damping 
term towards S(r) such that in the steady state, rj — > S(r). 
The advective derivative term (3 1 diTj ensures that the motion 
of the punctures, which are locally advected by /3\ is taken 
into account in the driving of r) —> S(r) (we want to drive rj to 
a specific behaviour in the neighbourhood of the punctures). 
To better interpret our suggestion for the dynamical gauge (|9j, 
it is useful to compare it with a simpler ordinary differential 
equation 



T Tt = - n + s{t) 



(10) 



Setting now r] = ij(t = 0), S a = S(t = 0), and assuming 
that T{dS/dt) <§; 1 so that it can be neglected at first order, 
Eq. ( |T0] > would have solution 



V - S+ (% - S )e 



-t/r 



(ID 



and thus ?; — » S as t — > oo. Although in the case of our gauge- 
evolution equation d9l, the source term is time-dependent and 
sometimes the time derivative can be rather large, especially 
near the punctures, equation ( fT0) i is useful to recognise that 
the damping term is itself damped and driven to the solution 
given by the source function S. 

The arbitrariness in the form of S(r) is removed in part by 
the works discussed in the previous sections and hence a first 
possible form is inspired by Q and thus given by 



Si(r) = J7mb(V) 



R ° (i-wy • (12) 



As we will discuss in the next sections, this choice works very 
well for single BHs, leading to smoother profiles for rj and 
consequently more stable evolutions. Similarly, another con- 
venient choice for the source is a combination of expression 
©and ((8) 



S 2 (r) = rjMB(r)r] S (r) 



Ro 



y/^diWdjW^ 

(i - wy , 



R 2 



r/s factor ensures minimal dynamics and implements the sug- 
gestion of Ref. 1 27 1 that rj — > at large radius to avoid the 
instability there. 

The generalisation of the source term (jT3j to the case of a 
binary system is straightforwardly given by 



S 3 (r) 



Ro 



y/fidiWdjW" 
(l-W) 2 



R 2 



R 2 



(14) 



where fj. , r*2 are the distances from f to the individual BHs 
centres. 

The following sections will be dedicated to the results of 
simulations performed with the evolution equation for the 
damping term d9b to study nonspinning single BHs via the 
source terms (fl2)-([T3|l, and to evolve BBHs via the source 
term ( |14) , As a final remark we note that although our dy- 
namical gauge requires the numerical solution of an additional 
equation [i.e., Eq. Q], the associated computational costs are 
minimal, given that we are evolving a scalar quantity and that 
the evolution equation is very similar to those already imple- 
mented within the Gamma-driver condition. More specifi- 
cally, the added computational costs range from 1 to 2% de- 
pending on the complexity and length of the simulation (the 
longer the simulation, the larger the impact of frequent restarts 
and thus the smaller the impact of the additional evolution 
equation). 



IV. APPLICATION OF THE NEW GAUGE: SINGLE 
BLACK HOLES 

We next examine the properties of the new dynamical 
gauge |9) for the damping factor and present its advantages 
when compared with the other position- and mass-dependent 
suggestions presented in the previous section. Specifically, we 
will study the t^mbM and rjs(r) prescribed forms for rj, and 
compare them with both the Si(r) and S^r) variants of the 
evolved gauge condition Q. We will start by considering the 
simple case of a single nonspinning puncture. 

Fig. [2] shows a comparison between an evolution using the 
prescription ?7mb (r) provided by Eq. (|3j (red dashed lines) and 
the new gauge using as source Si (r) from Eq. ( 12 1 (blue solid 



I \r 2 + R 2 J ' 
(13) 

This choice ensures that 77 is dynamically adapting in the in- 
ner region due to the tjmb factor, while in the outer region the 



lines). In both cases, the value of the damping factor at the 
puncture adapts automatically to the mass of the BH, due to 
its dependence on the conformal factor. However, when the 
damping term is given the form t/mb(»"), large spikes develop 
at the origin at late times and noise travels outwards as the 
gauge settles (cf. left panel of Fig.|2|i. These large noise pulses 
leave sharp features in r] which might lead to coordinate drifts 
and eventually affect the stability when simulating BBHs. On 
the other hand, when evolving the damping factor and using 
ilMs{r) as a source, it is possible to avoid any forcing of the 
damping term, which is instead always the solution of the dy- 
namical driver. This, in turn, reduces the noise and ensures 
long-term stability of the simulation as shown by the right 
panel of Fig. [2] 

We can further improve the stability properties of the damp- 
ing factor in the outer wave region by matching it with a 
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FIG. 2: Profiles along the z-axis of the damping term 77 for an isolated Schwarzschild BH. Different lines refer to the case when r\ is prescribed 
using Eq. l|5j (red dashed line) or when evolved in time using Eq. l|9j with source given by \ \ 2\ (blue solid line). The left panel refers to 
an earlier time and focuses on the central region of the grid. The right panel refers to a later time and shows a larger portion of the grid, 
highlighting that in contrast to the dynamical rj, the evolved one leads to smooth profiles at the patch boundary (see inset). Note that the outer 
boundary is at 400 M and cannot be responsible for the appearance of the spikes via reflection. 



function which drives it smoothly to zero. Fig. [3] shows a 
comparison between a BH evolution using the prescription 
VMB(i~)r]s(r) given by Eqs. |5]) and |8]l ( red dashed lines), 
and the evolved r\ using as source Eq. ( |13| > (blue solid lines). 
Note that in the case of the prescribed r/MB(r)r]s(r), the out- 
going gauge pulses still produces sharp features near the BH 
at z ~ 3 M, which are neither propagated away nor damped 
in-place. Moreover, sharp features far from the BH continue 
to be produced as the initial spikes pass through the interpatch 
boundaries. 

Smoother profiles can instead be obtained by evolving the 
damping factor, using rj MB (r)rjs(r) as a source. This is es- 
pecially true near the BH, while large variations but of small 
amplitude are still produced as the gauge pulses pass through 
the interpatch boundaries. Finally we note that in contrast to 
all the other evolution variables in our code, no artificial dis- 
sipation is imposed on the damping term so that the features 
of its evolution equation can be better appreciated. 



V. APPLICATION OF THE NEW GAUGE: BLACK-HOLE 
BINARIES 



As shown in the previous section, our evolved damping fac- 
tor leads to smoother profiles and consequently to stable BH 
evolutions, free of coordinate drifts. In this subsection, we 
study the effect of using different functional forms for the 
sources in the evolution equation Q. 



A. Equal-mass binaries 



We first consider the evolution of an equal-mass nonspin- 
ning BH binary, whose properties can be found in Table [H] 
For simplicity we have considered a system with small sep- 
aration D = 7 M, so that overall the binary performs only 
about 3 orbits before merging and settles to an isolated spin- 
ning BH after about 200 M. 

Fig. [4] shows a comparison of the profile of the damping 
term on the z-axis when using the evolution equation |9) and 
the source term given either by Si [cf. Eq. (12i; red dashed 



line], or by £3 [cf. Eq. ( 14 1; blue solid line]; in both cases we 
have set R = 20M. 

It is clear that when using either expression for the source 
function, the value of r\ at the location of the punctures adapts 
in time through the coupling with the conformal factor W, 
which tracks the position and the masses of the two BHs. It 
is also worth noting that near the two punctures, the two evo- 
lution equations yield very similar solutions for the damping 
factor as one would expect since R ^> D. However, when 



using the source ( 12 1 (red dashed line), the evolution of the 
damping term produces large gauge pulses which travel out- 
wards and are amplified by the mesh-refinement and the in- 
terpatch boundaries. This is particularly evident in the right 
panel of Fig. |4]which refers to a later time when the BHs have 
already merged. These undesirable features are very effec- 



tively removed by adopting the source term ( 13 1 (blue solid 
line), which provides a natural fall-off for the damping term 
as this propagates towards the outer boundary. 
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100 



FIG. 3: The same as in Fig. [2] but the prescribed damping term Eq. |5j is tapered with a function Eq. |8](, while the source term for the evolution 
equation is given by l |13| . Note that in this case the damping term falls off as expected but that in the case of the prescribed T]MB(r)r]s(r), 
the outgoing gauge pulses still produces sharp features near the BH at z ~ 3 M, which do not propagate away. Moreover, sharp features 
are produced as the initial spikes pass through the interpatch boundaries located at z ~ 40 M. Smoother profiles can instead be obtained by 
evolving the damping factor. 




FIG. 4: Profiles along the z-axis of the damping term r\ for an equal-mass BBH system. Different lines refer to the case when r\ evolved 
with sources given respectively by l |12| ( (red dashed line) or by \13\ (blue solid line). The left and right panel refers to two different times and 
highlight the importance of the fall-off term to avoid reflections at mesh boundaries. Note that the outer boundary is at ~ 2000 M and cannot 
be responsible for the appearance of the spikes via reflection. 



B. Unequal-mass binaries 

We next consider the evolution of an unequal-mass non- 
spinning BH binary with mass ratio q = 1/4, initial separation 
D = 8 M and which performs about 4 orbits before merging 



and settling to an isolated spinning BH after about 600 M. A 
complete list of the binary's properties can be found in Ta- 
ble [II] We find that an evolution equation for -q is a very con- 
venient choice also for unequal mass binary simulations. 
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Configuration 


mi/M m 2 /M 


xi/M x 2 /M 


Pf/M Pl/M 


single BH 


0.5 


0.0 — 


0.00000000 0.00000000 


BBH, q = 1 


0.5 0.5 


3.5 -3.5 


-0.00335831 0.12369380 


BBH, <7=l/4 


0.8 0.2 


1.6 -6.4 


-0.00104474 0.07293950 



TABLE II: Initial data parameters for the configurations studied. Expressed in units of the total mass M are: the initial irreducible masses of 
the BHs (obtained by iterating over bare mass parameters) rrii, the coordinates of the BHs on the a>axis Xi, the momentum of BH 1 PI (the 
momentum of BH 2 is equal and opposite, P 2 = —PI). The unit M is chosen such that each BH has mass 0.5Af in both the single and binary 
BH cases. 



1. Spatial Dependence 

Fig.[5]shows the spatial dependence of the damping term at 
some representative times. The two top panels refer to when 
the punctures are close to the s-axis for the cases where 77 is 
kept constant in space and time (red dashed line) and where 



it is evolved in time with a source given by ( 14 1 (blue solid 
line). Clearly, also in the unequal-mass case the new gauge 
drives the damping term to follow the BHs and to be smooth 
elsewhere. The bottom left and right panels show the damping 
term on the (x, y) plane at two times which are similar to those 
shown in the top panels and show the non-trivial but smooth 
distribution of the damping term which adapts to the motion 
of the punctures (the latter are near the maxima of the 77). The 
value of 77 is ~ 1.5/M near the large BH and ~ 4.5/M near 
the small BH, leading to a value for rrii rj ~ 1 for both BHs, 
where mi = 0.8 M and m 2 = 0.2 M are the irreducible 
masses of each BH. This demonstrates that the gauge con- 
dition is adapting the value of 77 to the mass of the BH. As an 
aside, we note that there is a region between the BHs where 77 
drops nearly to 0, and a "wake" of low 77 which follows behind 
the motion of the smaller BH. 



2. Waveforms and Errors 

One of the most useful quantities calculated from a BBH 
simulation is the gravitational waveform. The physical wave- 
form is that measured at future null infinity, and is indepen- 
dent of the gauge condition used for the evolution. In practice, 
waveforms are often computed at finite radii and extrapolated 
to future null infinity, and the gauge condition can have an 
effect on the extrapolation error. More importantly, different 
gauges can lead to different truncation errors in the evolution 
of the BH motion depending on how well we can reproduce a 
given field variable at the chosen resolution, and these errors 
will be reflected in the waveform errors. As a result, differ- 
ent gauges can in practice lead to small differences also in 
the calculation of the waveforms, and it is important to study 
the impact of the different gauges on the waveforms and their 
truncation errors. 

In Fig. [6] we show the real part of the I = 2, m = 2 mode 
of the gravitational waveform ^4 as extracted at r = 100 M. 
Different lines (dotted, dashed, solid) refer to the different res- 
olutions reported in Table [I] while the two panels refer to the 
cases where 77 is kept constant in space and time (upper panel) 



(lower panel). As expected, the differences in the waveforms 
due to the different truncation errors for the two gauges are 
very small and confirm that a different prescription for the 
gauges does not influence the gravitational-radiation signal. 
Figure [7] offers a different view of the gravitational -wave sig- 
nal by showing the phase difference in the I — 2, rn = 2 
mode, A022, of the waveforms as computed between the low 
and medium resolutions (red dashed line) and between the 
medium and high resolutions (blue solid line). The latter 
has been scaled to compensate for an 8th-order convergence 
which is indeed obtained as shown by the very good overlap 
between the dashed and solid curves during all of the inspiral, 
being slightly worse during the merger, when the convergence 
order drops. The left panel refers to the case when 77 is kept 
constant in space and time, while the right one refers to when 



and where it is evolved with a source given by S3 [cf. Eq. ( 14 1] 



it is evolved with a source given by ( 14 1 and R = 20 M. Once 
again the phase errors are comparable between the two gauge 
conditions. 

We see that adopting an evolution equation for the damping 
term allows one to reproduce with a comparable accuracy the 
numerical results obtained with the more standard prescrip- 
tion of a constant value for 77. At the same time, however, it 
also shows that in this way no special tuning is required and 
the prescription that works well for equal-mass binaries is also 
very effective for an unequal-mass case with q = 1/4. We ex- 
pect this to be true also for much smaller mass ratios, whose 
investigation goes beyond the scope of this paper but will be 
pursued in our future work. 



3. Impact on the Apparent-Horizon Size 

Some final consideration will now be given to the effect 
that the new gauge condition has on the apparent horizon 
(AH) coordinate size as it varies during the simulation. We 
recall that numerical simulations of moving punctures have 
routinely reported a certain dynamics in the coordinate size 
of the AH. This dynamics takes place during the early stages 
of the evolution, as the gauges evolve rapidly from their ini- 
tial conditions and reach the values they are brought to by the 
time-dependent drivers. Although these changes are of purely 
gauge nature and do not influence the subsequent evolution of 
gauge-invariant quantities, they represent nevertheless a com- 
putational nuisance as they require a high-resolution mesh- 
refinement box (i.e., the one containing the BH) to be suf- 
ficiently large so as to accommodate the AH as it grows. 
Clearly, this requirement becomes particularly important for 
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FIG. 5: Spatial dependence of the damping term r\ for an unequal-mass BBH system with q = 1/4. Different lines refer to the case when r\ is 
kept constant in space and time (red dashed line) or when it is evolved in time using Eq. l|9j with source given by (blue solid line). The 
top left and right panels refer to two different times when the punctures are close to the a>axis and show that the new gauge drives the damping 
term to follow the BHs and be smooth elsewhere. The bottom left and right panels, instead, show the damping term on the (a;, y) plane at two 
representative times and show the non trivial but smooth distribution of the damping term which adapts to the motion of the punctures (the 
latter are near the maxima of the r)). 



binaries with very small ratios, as in this case it is desirable 
to have the least dynamics and thus reduce the computational 
costs. 

To report the ability of the new gauge to reduce the varia- 
tions in the AH size, we show in the left panel of Fig. [8] the 
time evolution of the average radii of the unequal-mass binary 
on the (x, y) plane before and after the merger. We recall that 
during the inspiral and merger we follow three different AHs, 



two corresponding to the initial BHs (i.e., AHi and AH2) and 
a third one, which is produced at the merger and contains the 
first two (i.e., AH 3 ). We also note that the two initial AHs 
can still be followed for a certain amount of time after a sin- 
gle AH is found comprising the two. Shown in the left panel 
of Fig. |8]with red long-dashed (large BH), dot-dashed (small 
BH) and dashed lines (merged BH) are the radii obtained for 
the case in which 77 is kept constant in space and time. The 
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FIG. 6: Real part of the I = 2, m — 2 mode of the gravitational 
waveform \&4 for the unequal-mass black-hole binary. Different lines 
refer to different resolutions (see text for details) and the two panels 
refer to the case when r\ is kept constant in space and time (upper 
panel) or when it is evolved with source given by (lower panel). 



radii obtained when r] is evolved with source given by ( fl4] i 
are shown with blue solid (large BH), dotted (small BH) and 
short-dashed lines (merged BH). 

In each case, the individual AHs first grow in radius as they 
initially adapt to the chosen gauge, i.e., for t < 20 M. Af- 
ter this initial stage, however, the two gauges show a differ- 
ent behaviour and in the new gauge the AHs tend to main- 
tain their coordinate size more closely than in the constant-?? 
case (cf. blue solid and red long-dashed lines). This ability 
to conserve the original size is evident also after the forma- 
tion of the common AH at t ~ 520 M, where the new gauge 
shows an average radius for AH3 which is essentially constant 
for t ~ 200 M, while it grows of about 10% over the same 
timescale when using a constant 77 gauge. These effects pro- 
vide benefits for the mesh-refinement-treatment in our BBH 
simulations. Namely, they allow us to reduce the extent of 
the finest mesh-refinement levels containing the AHs and re- 
duce therefore the computational cost of the simulations while 
keeping a given resolution. 

Besides the changes in the overall coordinate size of the 
AHs discussed above, it is interesting to consider how much 
their shape changes in time in the two gauges and this is sum- 
marised in the right panel of Fig. [8] where we show the evo- 
lution of the ratio of the AHs' proper circumferences on the 
(at, z) and (x, y) planes, C xz /C xz , as computed for the three 
BHs (the first two panels from the top refer to the AHs of 
the binary, while the third one to the merged AH and thus 
a different time range). Overall, in both gauges the merging 
BHs remain spherical to a few parts per thousand. We have 
also verified that the masses of the BHs as computed from the 
AHs are consistent between the two gauges within numeri- 



cal errors. Interestingly, the new gauge leads to smaller os- 
cillations in the ratio for the smallest of the BHs. This is a 
very small improvement, which however minimises spurious 
gauge-dynamics and helps when making AH-based measure- 
ments. 



VI. CONCLUSIONS 

Even with a complete computational infrastructure, 
numerical-relativity simulations of inspiralling compact bi- 
naries would not be possible without suitable gauge condi- 
tions. A large bulk of work developed over the last decade 
has provided gauge conditions for the lapse and the shift 
which have been used with success both in vacuum and non- 
vacuum spacetimes when simulating binaries with compara- 
ble masses. However, as the need to investigate black-hole bi- 
naries with small mass ratios increases, evidence has emerged 
that the standard "Gamma-driver" shift condition requires a 
careful and non-trivial tuning of its parameters to ensure long- 
term stable evolutions of such binaries. 

As a result, a few different suggestions have been made re- 
cently in the literature to improve the Gamma-driver condi- 
tion for the shift and these have focused, in particular, on the 
specification of a spatially dependent damping term r\. This 
approach has been shown to work well under some conditions 
but not always and prescriptions which are effective in some 
cases can lead to instabilities in others. In addition, the pre- 
scriptions require the specification of coefficients whose tun- 
ing may be dependent on the mass ratio in a way which is not 
trivial. 

Following a different approach, we have presented a novel 
gauge condition in which the damping constant is promoted 
to be a dynamical variable and the solution of an evolution 
equation. We show that this choice removes the need for spe- 
cial tuning and provides a shift damping term which is free 
of instabilities for all of the spacetimes considered. Although 
rather trivial, our gauge condition has a number of advantages: 
i) it is very simple to implement numerically as it has the same 
structure of the other gauge conditions; ii) it adapts dynami- 
cally to the individual positions and masses of the BBH sys- 
tem and could therefore be used also for binaries with very 
small mass ratios; Hi) it reduces the variations in the coordi- 
nate size of the apparent horizon of the larger black hole thus 
limiting the computational costs; iv) all of the complexity in 
the new gauge is contained in the source function which can 
be easily improved further. This last point will be part of our 
future research in this direction. 
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FIG. 7: Phase difference in the £ — 2, m — 2 mode of the gravitational waveforms as computed between the low and medium resolutions 
(red dashed line) and between the medium and high resolutions (blue solid line). The latter has been scaled to compensate for an 8th-order 
convergence which is indeed obtained as shown by the good overlap between the curves. The left panel refers to the case when r\ is kept 
constant in space and time, while the right one to when it is evolved with source given by 
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FIG. 8: Left panel: Evolution of the average radius of the AHs on the (x, y) plane before and after the merger. Shown with red long-dashed 
(large BH), dot-dashed (small BH) and dashed lines (merged BH) are the radii obtained for the case in which 77 is kept constant in space and 
time. Shown instead with blue solid (large BH), dotted (small BH) and short-dashed lines (merged BH) are the radii obtained for the case in 
which 77 is evolved with source given by 1 14 1. Right panel: Evolution of the ratio of the AHs' proper circumferences on the (a?, z) and (a;, y) 
planes, C xy , C xz as computed for the three BHs. The first two panels from the top refer to the AHs of the binary, while the bottom one to refers 
to the merged AH and thus covers a different range in time. 
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